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Abstract 

We present a novel class of particle methods with deformable shapes that achieve high-order 
convergence rates in the uniform norm without requiring remappings, extended overlapping 
or vanishing moments for the particles. Unlike classical convergence analysis, our estimates 
do not rely on the use of a smoothing kernel but rather on the uniformly bounded overlapping 
of the particles supports and on the smoothness of the characteristic flow. In particular, 
they also apply to heterogeneous particle decompositions such as piecewise polynomial bases 
on unstructured meshes. In the first-order case which simply consists of pushing forward 
linearly transformed particles (LTP) along the flow, we provide an explicit scheme and establish 
rigorous estimates that demonstrate the convergence in the uniform norm and the uniform 
boundedness of the resulting particle overlapping. To illustrate the flexibility of the method 
we develop an adaptive multilevel version where particles are dynamically refined based on a 
local error indicator. Numerical studies allow to assess the convergence properties of this new 
particle scheme in both its uniform and adaptive versions, by comparing it with traditional 
fixed-shape particle methods with or without remappings. 

1 Introduction 

Efficient and simple particle methods are extremely popular for the numerical simulation of trans- 
port equations involved in many physical problems ranging from fluid dynamics [3 [T^] to kinetic 
(e.g., Vlasov) equations [111 US]. However, particle methods also suffer from weak convergence 
properties that lead to important disadvantages in many practical cases. Specifically, it is known 
that they only converge in a strong sense when the particles present an extended overlapping, that 
is, when the number of overlapping particles tends to infinity as the mesh size h of their initializa- 
tion grid tends to 0, see e.g., Refs. [31 [53]. Moreover, convergence rates are known to be suboptimal 
and to require additional constraints (namely, vanishing moments) for the particle shape functions 
that prevent high orders to be achieved with positive shapes. In practice, extended particle over- 
lapping can be expensive and it involves an additional parameter to be optimized, such as the 
overlapping exponent q < 1 for which the particles radius behaves like /i* ^ h. In Particle-In-Cell 
(PIC) codes for instance, taking q < I typically leads to increasing the number of particles per 
cell faster than the number of cells, since the latter determine the radius of the particles [12] ■ In 
Smoothed Particle Hydrodynamics (SPH) schemes it amounts to increasing the number of neigh- 
bors, i.e., interacting particles 23,. Because of these issues, many particle simulation methods do 
not meet the conditions of convergence which can result in numerically intensive simulations for 
acceptable results. Also, limited numerical resources often produce strong oscillations seen as a 
statistical noise that hamper interpretation of results and can further cause large scale errors. 

To suppress noise, many methods (like the Denavit redeposition scheme [15j . recently revisited 
as a Forward Semi-Lagrangian scheme (FSL), see e.g., Refs. [HKTOKIH]) use periodic remappings, 
i.e., particle re- initializations that smooth out the evolution. However, frequent remappings can 
introduce unwanted numerical diffusion which in many cases contradicts the benefit of using low- 
diffusion particle schemes, and to reduce this adverse effect some authors have introduced high- 
order adaptive remappings, see, e.g. Refs. [51 
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In this article, we present a new class of particle methods with deformable shapes that converge 
in the uniform norm without requiring remappings, extended overlapping or vanishing moments 
for the particles. Unlike classical error estimates based on a smoothing kernel argument, our 
analysis applies to general particle collections with Lipschitz smoothness and bounded overlapping 
properties. In particular our estimates easily extend to the transport of heterogeneous "particle" 
approximations such as standard finite elements bases. 

Our results are threefold. On a formal level first, we establish high-order convergence rates in 
L°° for a class of transport operators where the particle shapes are deformed with polynomial map- 
pings, the coefficients of which involve spatial derivatives of the backward flow. In particular, the 
first-order case is a linearly-transformed particle method (LTP) where the particles are deformed 
with the Jacobian matrices of the backward flow. It corresponds to a method already studied by 
Cohen and Perthame [5] who established its first-order convergence in but did not provide a 
numerical scheme for the deformation matrices. It can also be viewed as a modified version of 
Hou's method [18j where instead of using a global deformation mapping, each particle is trans- 
ported by the linearized flow around its trajectory. On a numerical level, we provide an explicit 
scheme for the LTP method that is based solely on pointwise evaluations of the forward flow, and 
we establish rigorous a priori estimates for both the transport error and the particle overlapping. 
To illustrate the flexibility of our approach we also present an adaptive multilevel LTP scheme 
where local estimates for the single-particle transport errors are used to decide which particles are 
dynamically refined, and a local correction filter for high-order positivity-preserving hierarchical 
approximations is presented. On a practical level we eventually present numerical results obtained 
with academic test cases. They show that a B-spline based LTP scheme converges with better rates 
compared to the traditional smoothed particle method (TSP) with extended overlapping. Here the 
convergence is sometimes improved by introducing periodic remappings, but compared with the 
FSL scheme we also show that optimal remapping frequencies are much lower with LTP, leading to 
lower numerical diffusion and computational cost. Finally we verify that our adaptive LTP scheme 
enhances the convergence of solutions with sharp edges by equi-distributing the transport error. 

As we have pointed out, deforming the particles is not a new idea. Dynamic adaptation, 
i.e., refinement and coarsening of particles, also has a substantial history. In vortex methods for 
instance, Cottct, Koumoutsakos and others have introduced a variety of algorithms to handle 
particles with spatially varying sizes based on global or local mappings, see e.g. [TT1|4]. In kinetic 
(PIC) schemes, Lapenta, Assous and others developed "rezoning" algorithms to increase or lower 
the number of particles per grid cell, while preserving the corresponding grid moments such as 
charge, current or energy density, see e.g. |201 [Tj. And for pure transport problems, Bergdorf and 
Koumoutsakos [5J have studied a wavelet-based FSL scheme with adaptive, high-order remappings. 
However, although the list is not exhaustive we observe that these methods adapt the size of 
particles with fixed shape, and hence require high frequency remappings. A noteworthy exception 
is the Complex Particle Kinetic method developed by Bateson and Hewett for the simulation of 
plasmas [51 [T7] where in addition to having their Gaussian shape transformed by the local shearing 
of the flow, particles can be fragmented to probe for emerging features, and merged where fine 
particles are no longer needed. When mature, our adaptive scheme should be compared with the 
above methods. 

The outline of the article is as follows. In Section [2] we begin with a rapid overview of the main 
particle methods, introduce some notations and state our main results. In Section [3] we present a 
class of high-order particle methods with polynomial deformations, as well as a fully discrete LTP 
scheme, and we establish a priori estimates. Section [4] is then devoted to the description of an 
adaptive multilevel LTP scheme based on refinable B-splines, and numerical results are presented 
in Section m 

Readers mostly interested in the practical aspects of the LTP scheme may prefer to first read 
Sections |2.3[ |3.1[ |3.3| and |3.4[ Its main properties are stated in Theorem |3.10[ and its adaptive 



multilevel version is summarized in Algorithms 4.1 and 4.6 
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2 A brief review of particle methods 



To introduce sonic notations and state our main results we begin with a rapid overview of particle 
methods. Following [9] we consider the linear d-dimensional transport equation 

dtf{t,x) + u{t,x)-Vf(t,x) = 0, ie[0,r], x eR'^ (2.1) 

associated with an initial data /'^ : M"^ M, a final time r and a velocity field u : [0, r] x ^ M"^. 
In fluid problems for instance, we have d = 2,3, while in kinetic formulations M.'^ is a phase space 
with d < 6, and u is a generalized velocity field with components of velocity and acceleration. We 
assume that u is smooth enough for the characteristic trajectories X{t) = X{t, tQ, xq), solutions to 

X'it)=u{t,X{t)), Xito) = xo, (2.2) 

to be defined on [0, r] for all xq E M'' and to G [0, t], see, e.g., [23]. In particular, the corresponding 
characteristic fiow Ft^^t ■ xq m- X{t) is invertible and satisfies Ft^to = {Fto.t)~^- Solutions to (2.1 1 
read then 

fit,x)^ f{to,{Ft,,ty\x)) foraU to, t e [0, r] and a: e M'^. 

For simplicity, we restrict ourselves to the incompressible case where divu = 0. In this case, the 
flow is measure preserving in the sense that its Jacobian matrix Jp^^ ^{x) — (^dj{Ftg t)i)i^^ has 
a constant determinant equal to 1, 

det {Jpt^^t {x)) = 1 for ah to, t S [0, r] and x G W^. 



2.1 The traditional smoothed particle method (TSP) 

In the standard "academic" particle method [53], numerical solutions are typically computed as 
follows: considering deterministic initializations for simplicity, the initial data f'^ is first approxi- 
mated by a collection of particles on a regular (say, cartesian) grid of step h > 0, 

fhA^) E Mn^eix ~ Xk) with xl hk 

and with weights 

Wk[f):^ I j\x)Ax or w^if) := f{xl). (2.3) 

Here = f-~'^f{'/f) is a particle shape function with radius proportional to e, usually seen as a 
smooth approximation of the Dirac measure obtained by scaling a compactly supported "cut-off" 
function ip for which a common choice is a B-spline. Particle centers are then updated at each 
time step t" = nAt by following the fiow = Ff,i in+i or its numerical approximation, and the 
weights are kept constant, leading to 

E MfWeix - xl+^) « f{e+\x) with xl+' := F-{xl). 

In the classical error analysis [3] [53], the above process is seen as (i) an approximation (in the 
distribution sense) of the initial data by a collection of weighted Dirac measures, (ii) the exact 
transport of the Dirac particles along the fiow, and (iii) the smoothing of the resulting distribution 
'^k{f^)Sx^ with the convolution kernel ip^. The classical error estimate reads then as follows: 
if for some prescribed integers m > and r > 0, the cut-off if has m-th order smoothness and 
satisfies a moment condition of order r, namely if J ip ~ 1, J\y\^\'p{y)\ dy < oo and 

J yl' . . . y"/ip{yi, ...,yd)dy^O for s e N"^ with 1 < si + • • • + Sd < r - 1, 

then there exists a constant C independent of /i or e, such that for all n < r/At we have 

ii/r) " /LIU- < c^e^wfww-^ + {h/er\\n\w--^) (2.4) 
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with 1 < n < oo. More recently, Cohen and Perthame [5] observed that defining the weights as 

Mf) ■■= I f{x)Mx-xl)dx (2.5) 



with a weighting function (ph = (pi-jh) derived from a continuous and compactly supported (p such 
that 

^ fcj^ . . . k^^ipiy -k) ^y{^ ... y^/ for s e N"* with < si H h Sd<m-1, 

one has the improved estimate 

Wfin - /LIU- < C[e^f\\w^.. + {h/erWfh.) (2.6) 



with a new constant that is again independent of /°, h or e. Note that (2.6) is better than (2.4 1 
in that m is not constrained by the smoothness of which allows to reach higher convergence 
rates. Indeed balancing the error terms in the above estimates suggests to take e ^ /i' with 
q = yielding a convergence in h'' = h . In particular, if /° G W'^''^ for some integer s 

then the best possible rate with standard weights is only obtained with m = r — s. 

With the improved weights instead, one can take a higher value for m and obtain estimates close 
to w-f"- Moreover, the latter approach also allows to improve (i.e., reduce) the particle 

overlapping, since the corresponding exponents are 9=5 and ~ 1, respectively. In either 
case, we see from the term h'^^\\f'^\\w^-'^ in the estimates that extended particle overlapping does 
not only make the simulations more expensive, it also deteriorates their convergence order. 

2.2 The forward semi-lagrangian scheme (FSL) 

In Forward Semi-Lagrangian schemes (the new name for the periodically remapped particle method 
introduced by Denavit [T5]). extended overlapping is usually not required and particles have the 
same scale than their initialization grid, i.e., e — h. Instead, their weights and centers are re- 
initialized using a regular grid once every time steps. Thus, denoting by 

Ah-, gi-^ ^ Wk{g)(ph{x ~ x°) 



the particle approximation operator with again x° = hk and weights computed, e.g., as in (2.3), 
and by 

r":<^^(._4^)^^,(._F"(a;^)) 
the fixed-shape particle transport operator, the FSL scheme takes the form 

/r^= E-rV.(---r^):=r"/^ with ~r,:=[y'^ -^^n^N^n ^2.7) 

\ fh otherwise, 

and where has been extended to collections of particles by linearity. 
2.3 The hnearly-transformed particle method (LTP) 

In this article we shall develop a lesser-known approach already studied by Cohen and Perthame [3] 
who observed that by transporting the particles with the linearized flow around their trajectories, 
one obtains a convergent method (in L^) with particles scaled with their initialization grid, and 
no remappings. On a formal level this amounts to defining linearly-transformed particles as 

pn,k{t,x):=p,,{Mt){x-Xk{t))) with |?S!J (2.8) 



In practice, we observe that occasional remappings are needed for accurate solutions. We may then 



rewrite the numerical LTP scheme in a form similar to (2.7), but particles are now associated with 
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invertible d x d deformation matrices representing backward Jacobian matrices at a;". Thus, 



numerical solutions read 



(2.9) 



and transporting the particles consists in updating the deformation matrices together with 
the particle centers x^, initialized as := Id and :— hk, respectively. In Section 3.3 wc will 



describe a numerical method for doing so, that is solely based on pointwise evaluations of the 
(approximated) forward flow F". We may summarize our findings as follows. 



Main results. The formal LTP method (2.8) converges with order 1 in L°° , and arbitrary orders 



can be reached with proper polynomial deformations which coefficients involve the derivatives of 



the backward flow (see Theorem 3.2). On the numerical side, an explicit implementation of the 
LTP transport operator based on finite- difference approximations of the forward Jacobian is also 
shown to converge with order 1 in L°° with no remappings required (see Theorem 3.1(J^ . In practice 
this leads to improved convergence compared to both TSP and FSL schemes, with lower remapping 



frequencies than the latter (see Section 5.1). Moreover, local error estimates for the single-particle 
transport errors are established (in Theorem 3.1(^ that can be used for dynamic refinements in a 
multilevel particle framework (see Section 5.2). 



In the sequel it will be convenient to use the maximum norm |lx||oo maxi|a;i| for vectors and 
the associated m||oo ■— max^ ^I^jl^ij | for matrices. For functions in Sobolev spaces W™-'°°{ijj) 
with Lo C W'' and integer index m > 0, we will use the semi-norm 



a a 



(2.10) 



and for conciseness we drop the domain when it is the whole space. 



3 Particle methods without smoothing 

In this section we present a class of particle methods that deviate from the "smoothed" approaches 



described in Sections |2.1| and 2.2 in the following sense 



• Convergence (including high-order) is proved without resorting to a smoothing kernel argu- 
ment as with TSP schemes. Instead particles have their radius proportional to the meshsize 
h of their initialization grid, which 

(i) permits using not only one single particle shape but also heterogeneous collections such 
as standard finite element bases, as no vanishing moments or high-order smoothness is 
required ; 

(ii) allows to compute their weights with standard approximation schemes ; 

(iii) yields uniformly bounded numbers of overlapping particles. 

• Remappings are not required for convergence as in FSL schemes, although they may improve 
the results in practice. 

To simplify the presentation we focus on homogeneous spline particles, although heterogeneous 
bases could be used, see Remark |3.3| below. Thus, in Section |3.1| we first review one local ap- 
proximation scheme with B-splines. In Section 3.2 we introduce particle transport operators T^^) 



with polynomial shape transformations that converge without smoothing arguments, and in Sec- 
tion 3.3 we describe one explicit implementation of the first-order operator T^^) . In Section|3.4 |we 
give a practical tool for localizing particles with linearly transformed supports, and in Section |3.5| 
we establish rigorous error estimates that prove both the uniform convergence and the uniformly 
bounded overlapping of the deformed particles. 
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3.1 High-order quasi-interpolations with B-sphne particles 

For simplicity, we consider homogeneous particles based on tensorizcd B-splines. Specifically, the 
univariate cardinal B-splines can be defined recursively as 

, 1 

'^o:=Xr 1 1, and Bp{x) := {Bp-i * Bq){x) = / ^ Bp-i for p > 1, 

where x denotes the set characteristic function. Clearly, Bp is a piecewise polynomial of degree 
p, it is globally C^^^ and supported on [— Moreover its integer translations span the 
space of cardinal splines of degree p, see, e.g., [H]. Thus Bi{x) — max{l — |a;|, 0} is the traditional 
"hat-function" , B3 is the centered cubic B-spline, and so on. Our reference particle shape function 
will be the tensorized, centered B-spline of odd degree p, 

d 

^{^) '-^W^pi^i) supported on supp(iy9) — [— Cp, Cp]'' with Cp := (3.1) 

i=l 

Note that B-splincs of odd degree are refinable on dyadic grids, which will be used in Section [1] to 
reformulate our particle method in a multilevel framework. 

Since our particles are scaled with their initialization (or remapping) grid, we can use standard 
approximation schemes that rely on the fact that the span of their integer translates 

^lA^)^^h{x^xl) = h-''^{h-^x-k), k&Z", (3.2) 

contains the space Vp of polynomials with coordinate degree less or equal to p. Specifically, we 
shall consider quasi-interpolation schemes described by and 25J, where high-order B-spline ap- 
proximations are locally obtained by pointwise evaluations of the target function. In the univariate 
case they take the form 

"^h"^^ ■ 9 '~^^Wk{g)ipl^k with normalized weights Wk(g) := h ^ a; 5(2;"+;) 

fceZ l'l<mp 

and symmetric coefficients a; = a^i defined in such a way that A^^'^^ reproducts the space Vp^'^\ 
They can be computed with the iterative algorithm in [51 Section 6] : for the first orders we obtain 

• TOp = and ao — I for p — 1, 

• nip — 1 and (oq, ai) = (|, — g) for p = 3, 

• nip^ A and (00,01,02,03,04) = (fifi "Hin' 2k' 3M0' THoo) P = 5. 

In the multivariate case we can tensorize the above, as it is easily checked that the operator 

Ah ■■ g ^ ^ Wk{g)'fl^k with Wk{g):=h'^ ^ o;g(a;2+,), 0;:= J| a;, (3.3) 

keZ'i l|i||=o<mp l<i<d 

reproducts any polynomial tt e Vp. Moreover, we have 

\\Ahg\\L^ < {2cpfy,,\\L^ snp\wk{g)\ < (2cp)''|l^|U^ |lo|l,i |l.g||ic. (3.4) 

with ||a||£i = <m l*^'!) by using the fact that no more than (2cp)'' B-splines overlap. It 

follows that Ah is uniformly bounded in L°° , with 



\Ah\\ciL-^ :=sup "f^f < (2cp)'*||^|U=e||a||,i. 



Using a localized version of (3.4), we write that for an arbitrary n £ Vq, q < p, 

< {\\Ah\\c{L°-) + l)ll3 " 7r||L=(B„(2;0,h(mp+Cp))) 

where Bca{x,p) denotes the open cube of center x and radius p. Taking for n the q-th Taylor 
expansion of g around Xk , we thus find for all q < p, 

\\A,,g~g\\L^ <h'i+^CA\gU+i with ca = {\\Au\\ciL^) + l Y ^^^^^j, . (3.5) 
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3.2 High-order particle transport with polynomial shape transforma- 
tions 

We now address the problem of transporting a collection of particles 

along the flow F — Fq^t associated with the whole time domain, with no remappings. Here we 
consider particles initially centered on the cartesian grid (3.2), with weights satisfying 

\wk\<c^h''\\f\\L^, keZ" (3.6) 

with a fixed Cw > 0, as it should be with any standard approximation scheme (if /° = A^f'^ we can 
take Cw = \\a\\gi). Moreover, to simplify the analysis we begin by forgetting the time discretization 
and consider that F is known and can applied exactly. (In Section 3.3 we will take into account the 
time approximation errors to study a discrete particle scheme.) Thus, in the traditional method 
particles keep their shape and are simply translated with 

T = 7(0) : (fl f, Lph{- - Xk) where Xk := F{xl). (3.7) 

Since the exact transport operator reads 

T,^:9^9{F-\-)), (3.8) 



the approximation (3.7) is exact for point (Dirac) particles. For finite-size particles however, the 
method does not converge. Take indeed p = I (i.e., (p is the standard hat function), and consider the 
infinitely smooth problem where /"^ = 1 and u{t,x) = {—X2,xi) in two dimensions, over the time 
interval [0, j]. Then any reasonable initialization will give Wk = h? , hence, f^{x) — 1, and clearly 
the exact final solution is /(f , x) = 1. Now, at the final time the particle centers will have rotated 
of T = ^, therefore every particle with |fci| + \k2\ = 1 will be centered on (cos(6' + ^),sin(6' + j)) 
with 9 £ "IN, and hence contributes to a; = with T(o)(p5J ^,(0) = h~^{l ~ ^)^7 in addition to iph,o 
which does not move. Since the other particles do not contribute to x = 0, the final error satisfies 



ll(T(o) - > 1^(0) A°(0) - 1| = 2(72-1)2, regardless of h. 

To improve the accuracy of the transport operator we shall estimate its error, using as basic 
property that the transported particles are localized by a relation 

suppiTifilj.) C Boo{xk,hph,k), (3.9) 

with radius factors uniformly bounded by p := sup^j^Q.feeZ'' Ph,k < oo. This clearly holds for T(q), 
taking ph^k = Cp = p. If T is such that (3.9) holds, then the particles overlap in a bounded way. 
Indeed for aU x,k such that Tip^ ^{x) ^ we have - h^^F^^{x)\\oc < /i~"^|-F~"^|i||xfe — a;||oo < 
Ph,k\-^~^\i' hence we may define an overlapping constant that satisfies 

e := sup sup #({fc e : T^lki^) ^ 0}) < {2p\F-'Uf. (3.10) 

h>Q xeR'' 



Moreover, it is easily seen that the particles transported with the exact operator (3.8) also have a 
bounded overlapping constant: if fc is such that Tex^/J^ ^^(a;) ^ then obviously \\ k— fi^F^^ {x)\\oci < 



Cp, hence 



Oex := sup sup #({fc e : T,^p>\^{x) ^ 0}) < (2cp)" 



In particular, we can decompose the global error (T — Tcx)fjl in terms of local (single-particle) 
transport errors eh,k = [T — Tcx)^Ph fc, as 

||(T-rex)/,°||L- <esup|lzi'fce;,,fc|lL- </j''c„|l/°|lL-esup|le,,,fc|lLoo (3.11) 
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holds with := + Qex, where we have used (3.6). We next locaHze the particles transported 
with the exact operator by writing 

T,^V^l,{x) ^ =^ F-i(x) e := supp((pO_fc) = B^ixlhcp) 

=^ \\x-Xk\\oo< l^^li,so JI^"^(a;)-a;°lloo < /iCp|F|i_EO^, 



so that 



supp(Te,^° ) = FiJ:lf,)^F{B^{xlhcp)) C B^ixk,hcp\F\,^^o ). 



(3.12) 



Hence the localization supp(eft_fe) C Bao{xkThph^k) for the error terms, with a new factor pfi^ '■= 
max{ph_fc, Cp|F|]^ 2j ^} ^ P '■= niax{p, Cp|F|i}. In the case of particle translations, i.e., for T = TJ-q), 
it follows that 



|e/i,fc||Loo = sn]y^^B^^s,^,h~pu.k)\Vh{x - Xk) ~ iph{F-^{x) ~ xX)\ 

< |<y9ft|iSup^gB^(5^ ,,^^^)||(F-i(a;) -.t) - {F-^{xk) - Xk)\\c 

h-'^-'^\(p\i, see (lO). Thus, we find 



(3.13) 



where we have used the scaling \(pii 



||(r(o) - Te.)/°|U^ < c^l^liOplF-i - /|i||/°|U~ 

for the fixed-shape particle method. Note that we can take Q = (2cp)''(l + |i^~^|f) and p = Cp\F\i 
in the above estimate, indeed det(Ji?(a;)) = 1 valid for all x gives 1 < || JF(a;)||oo < l^^li- Not 
surprisingly, the above analysis fails short of proving the convergence of the method. But it 
suggests one way to improve its accuracy: letting indeed 



0fc(s) = (t)k{s; x) (F ^ - I){xk + s{x - Xk)) 



(3.14) 



we may see the approximation {F ^{x) — x) ^ {F ^{xk) — Xk) involved in p. 13 ) as a zero-order 
Taylor expansion (/'fe(l) ~ 0^(0), and consider replacing ipii{x — Xk) by 'Ph{^k.(r){x)) , r > 1, where 



$,.,„(x) ■.^x-Xk + MO) + --- + -/k\0) ~ *fe,ox(a;) := F'^x) - xl 



(3.15) 



is formally an r-th order approximation. One could also consider expansions of the alternate 
(^fc(s) := (J - F){xl + s{F-'^{x) - xD) since 0^(1) - 0fc(O) = (t)k{l) - 0fe(O), but the particular 
form of (3.14) gives 



c^^k\^) = T.---T.h---9iAF-'-I)i 



Xk + s(x - 



^k))l[{ 



X - Xk)u 



i=l 



(3.16) 



so that ^k,{r) is a polynomial mapping which coefficients involve evaluations of space derivatives 
of F^^ at Xfc, which can be written in terms of the derivatives of F, at F~'^{xk) = x^. Also, (3.16) 
allows to specify the accuracy of the Taylor expansions (3.15) as follows: for all r > 1 and every x 
in a localized domain u C Boo{xk, hp) with p > 0, we have 



where (lu) denotes the convex hull of u 



< h 



r+l (P) 



r+1 



(r + l)! 



(3.17) 



For r = 1 we observe that Jp-i (xk) — (Jpixl))^^ is the Jacobian matrix of the backward flow, 
or equivalently, the matrix inverse of the forward Jacobian. The linearly-transformed particles are 
then defined as 



T{i)fl,k{x) := V'h{^k,(i){x)) = tphiJkix - Xk)) with 



Xk := F{xl) 
Jk-.^ {Jf{xI))-\ 



(3.18) 
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Since <&fc.(, ) is invertible the particles have a parallelogram support localized with 
supp(r(i)(^^ J.) = Xk + Jk^iBoo{0,hcp)) C Boo{ph,k,hcp\\Jk^\\oo), 



i.e., (3.9) now holds with ph^k — Cp\\{Jk) ^||oo < — p. As a consequence the above analysis 

applies readily to T(^i), the only noticeable change being that instead of ( 3.13^ we now write 



< |¥'h|iSup^gB^(j^ ;,^^_^)||$fc,(i)(a;) - $fe^ex(a;)||oo 



(3.19) 



with pii^k = Cpmax{|l(Jfc) ^||oo, l-^li.s° ^} — Cpl-F'li.E" Using (3.11) and (3.10) we thus find 
11(^(1) - Te,)/0|U^ < /if ec^|^|i|F-i|2||/"|U^ 

for the formal LTP scheme, with p = Cp\F\i and 8 = (2cp)''(l + \F-^\f\F\f). 

Remark 3.1. When applied to a particle with general shape and support containing x^, the LTP 
transport operator is defined as the exact transport for the linearized flow 

F,o :x^Fixl) + JFixl){x-xl), 

i.e., we set T^D'flk - T,^[Fx<l]^U = ^h,k ° iP^ir'- 

For larger values of r, care must be taken when defining T(^r)- Indeed, ^k,{r) may not be 
invertible hence there is no guarantee that \\x — Xk\\oo ^ h in the support of (y9/i(<i>j, (r)). To 
overcome this we define the transported particle as being restricted to some domain Jlh.k that is a 
priori bounded, i.e., we set 



T{r)vl,kix) = iphi'i>k,(r){x))xt^Jx), 



(3.20) 



where again ^k,{r) is the polynomial mapping defined by ( 3.15 )-( 3.16) and x denotes the set 
characteristic function. To determine S/i.fc we next observe that in order to carry out the error 
analysis, a convenient property is that 



Fi^lk) C S,,,fe C B^{xk,hph,k) 



for some ph,k- Indeed, from (3.12) this would give supp(eh,fc) C ^h,k and 

\\eh,k\\L-^ "^^^Pxeth.k \'Ph{'^k,{r){x)) - fh{'^k,cx{x))\ 
I'/'ll ' (r+1)! K '\r+l,B^(xk,hph,k)- 



Using (3.12) we see that an obvious solution to (3.21) is obtained with 

'^h.k ■= Boo{xk,hphM) with ph,k := Cp\F\^^^o^ ^ 



(3.21) 



(3.22) 



(3.23) 



However, one may want to restrict the particles to parallelogram supports that are closer to J^(S° j.) 
than the upper bound Bao{xk,hph,k)- To this end we observe that for any x G F(Y1'^ ^), applying 
(3.17) gives 



||*/c,(i)(a;)||oo < \\F ^{x) - xlWoo + \\^kAi){x) - $fe,cx(a;)|loo < hcp{l + A) 



as long as 



A> ^hcp\F\l^^Q^jF ^|2XF(so_j)- 



It follows that a finer solution to ( 3.21[ ) is given by 

(Eh,k ■■= ($fe,(i))"'(Soo(0,/iCp(l + A))) =ife + (Jfc)-i(Soo(0,/iCp(l + A))) 
\ph,k ■■= Cp|F|i ^0^(1 + A). 

In either case, global estimates are easily derived from the above discussion. 



(3.24) 
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Theorem 3.2. As above we consider B-spline particles (3.1|-(3.2| initially centered on regular 
nodes, and denote by F = i^o.r the exact characteristic flow associated with the transport equa- 



tion (2.1 1. We let Tf^r) be the r-th order particle transport operator defined as (3.18) for r = 1, or 
as (3.201 for r > 2, with restriction domain given by (3.23) or (3.24|. In the latter case we take 
A > 2^Cp|-F|^|-F~^|2, otherwise we set A := 0. Then the transported particles are localized by 

supp(T(,)(pO_J C B^{F{xl), hp) with p = Cp\F\i{l + A) 

and their overlapping constant satisfies 

e sup #{{k e : T(r)AM ^ 0}) < (2p|^^^'|i)'- 

Moreover, the global transport error satisfies 

\mr) Te.)/°||i^ < h^c^\vUe^l^\F-\+,\\f"U^ 



where G) ~ (2cp)'* + O, and where c^ is the constant from (3.6) 



Remark 3.3 (heterogeneous "particle" approximations). As previously pointed out, our arguments 
do not rely on a smoothing kernel unlike classical analysis of particle methods [21 [53], and the above 
estimates readily extend to the heterogeneous case where the "particles" iph,k are not derived from 
a reference Lp but instead are defined as piecewise polynomials with global continuity constraints 
(e.g., standard finite element bases) on unstructured meshes of M"^, under the usual shape regularity 
and quasi-uniformity assumptions. 



3.3 A finite diff'erence implementation of the LTP method 

We now describe and analyze an implementation of the particle transport operator T^i-^ that only 
requires pointwise evaluations of some forward numerical flow 



Fit = K„ 



(3.25) 



given by an explicit solver for the ODE (2.2) over the time step Specifically, we consider 

the following approximations. First, the exact transport operator 



rpn rp r pn ] 



(3.26) 



when applied to a linearly transformed particle "ysJJ j, = iph{D'^{- — x^)) as in (2.91, is approached 
by the LTP transport operator according to Remark |3.1[ i.e.. 



T^,)[F:^J : iplk ^ ^h{DIJI,^[- ~ F^^ixl))) with Jl^^ ^ {J f^^M))'' ■ 

The Jacobian matrices involved in the latter are then approached by a finite difference scheme, 
and finally the values of F^^ are replaced by those of the numerical flow F". Thus, we set 



h 



)) 



(3.27) 



with x"^^ := F"-{x^J.) and -DJ!^^ D^Jk- Here w {•^F;^^{xk))~^ is a numerical approximation of 
the backward Jacobian matrix obtained by first approximating the forward Jacobian by a centered 
finite difference scheme, 



he-i 



{F-Uxl - he,)) ^ d,{F-^Uxl) 



(JD.,, (2/.)-i((F"),(.xI! 
for i, j — 1, . . . , d, and then letting 

:= det(J^')3(J,")-i or simply J," := {Jj:)-\ 



(3.28) 



(3.29) 



whether we want a conservative transport (i.e., where J TJ^^^ k ~ I fc)' '-'^ 

Note that det{Jp^^) = 1 on M'^, therefore it is reasonable to assume that the dx d matrix JJ} is 
invertible. In the following Lemma we establish a sufficient condition for this, together with some 
a priori estimates for the resulting approximations. 
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Lemma 3.4. Let e^^ := — -F^Hl^ be the error of the ODE solver, and write 

where B^f^^^ := Booix^,h). Then the approximated determinant is bounded as 

det(Jfe") > 1 - dV2,JMM)'"'- (3.30) 

Moreover, in both the conservative and the non-conservative cases, the approximated Jacobian 
matrix JJ} defined in (3.29) satisfies 

\\Jk - V."j-(^c" (^fc))lloo < 2d'h^,l,{9'^r^^^l,)'^''-^^ 

and 

provided 0j! > 0. In particular, for all p > we have 



sup 



(3.31) 
(3.32) 

(3.33) 



with 



Proof. For conciseness, we denote 



ph)- 



2h 



l<i,j<d 



The latter is the finite difference approximation of J^'™ obtained by substituting F" by in 
(3.28), so that \\Jj!'* — JkWoo < h~^de'^i. We observe that with the semi-norms (2.10) we have 



,,\\Jk'*\\oo < l^cxli^s;:, hence 



^^fc 1 1 oo < 1 1 t/fe 1 1 oo 



jn _ j7i 



< \F:i,\,^Bi,+h-'del,^pl, 



(3.34) 



Next we write two Taylor formulas for s i— >■ F^^{x^^ + scj) with j = 1, . . . ,d, namely 



ah 



F^Jxl + ahe,) = F^^{xl) + crhd^F^^xl) + {ah- s)djF:^^{xl + se,) As, a = ±1 

JO 

By taking their difference we obtain \\J^'* ~ >/fe'™||oo < ^^l-FScb.s;; , hence 



7" _ 7" 



<K\\Fll\2,B^^^+h-'de\,)^ pl 



(3.35) 



Using det( J^''^'^) = 1 and the multilinearity of the determinant, we then find 

|det(jn - 1| < d\\J^ J^''=''|Umax{||J,"||^,||J^'^''|U}''-^ < dhpl,{pl,f-' (3.36) 



which shows (3.301. From the cofactors formula A ^ — det(yl) ^C* involving the transposed 
cofactor matrix we next write for 9^ > Q 



det(J^ ) 



where we have used J^'°^ — iJf^'°^) ^- Similarly, we have 

II n'-wd-i 

\\{Jk)-'\U < d^-j^jf- < d{eir\pi^Y-^ 

det(Jfc j 



(3.37) 
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hence in the non-conservative case, the approximated Jacobian satisfies 



'k l|oo||t^fe 

For the conservative case, we infer from 6*]! < min{l, det(J^')} that for a < 1 

^dot(J^) 



|l-det(jn"| 
Taking a — ^, we then obtain 



az 



dz 



<H(0^)"-^|l-det(Jfc")|. 



(3.38) 



, - jr'^iioo = iidet(j;')3(j,")-i - {Jin-'w^ 

< idet(jfe)3 - i|ii(j^)-i|i^ + ii(j,")-i|uii J." - jr'^iuii jr'^iic 

<dh^^l,{^^l^f'^'-'\{elY^-^ + d{el)-^). 



Estimate (3.31) is then easily derived from the above inequahties, using d>l and < 1. Next 
we note that (3.35) gives ||(J^)~^ — {Jk''^^)~'^\\oo < h^^k non-conservative case, while in 

the conservative case it gives 

11(4")-' - (4""'=)-'lloo = ||det(J,")-^ 4" - jm\oo 

<|det(Jn-^-l||lA"lloo + ||Jr-Jr''l|oo 
<hlilkil + {9l)-'^-\t^lkY) 



where we have used (3.381 with a = — together with (3.36), (3.34). This shows that (3.32) is 
valid in both cases. To complete the proof we finally write 

+ 1! j^-°^(F"(x;') - F:^{xm\oo + II - m^D) - ((^c" )-'(^) - ^fc)iioo 



where we have used 
one. 



•^fc^' ~ ^"'(.^k) first inequality and (3.17) with r = 1 in the second 



□ 



3.4 Intermission: localization of particles with deformed shape 

In order to compute pointwise values of the density fj^{x) — ""^fc fe(^)j it is important to 
access every particle which support contain any given x £ M'* in a reasonable amount of time. This 
can be done with a pre-processing algorithm that subdivides the phase space into simple domains 
such as dyadic cells fi™ = 2~^{m + [0, 1)''), m e Z'^, with 2~^ sa h, and then writes in the scope of 
every such cell the indices of the overlapping particles, namely 

Kl^ikeZ" ■.^l„(n^)^0}, meZ"^. (3.39) 

In practice, one can run a cell marking algorithm for each unstructured particle (^^ , starting with 
the cell fim containing x^, i.e. m — [2-'a;j!j, and recursively testing whether the adjacent cells 
overlap the parallelogram support of <pJJ j, . To perform this test we can use the following result 
with A = {D^)-\ B = Boo{0,hcp) and S' = Boo(2-J(m + i) - x^, 2-J-i). 

Lemma 3.5. Let B — HiLiI^^i ~ TiTXi + r^] and B' = JliLiI^i ~ ^iJ^i + ''il ^'^^ orthotopes 
aligned with the coordinate axes, and A an invertible d x d matrix. We have 

ABnB'^tll <^ liAx - x'),\ < fi + r'^ and \{x - A~^x')^\ < r, + for i = l...d 
with h ■■= Ej=il^i,jki and f[ := E^=ilA\i kj - 
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Proof. Observe that S := l\'^^^[(Ax),-n,{Ax)i+n] andS' := Il^^i [(^"^a;)i-fj, (A-ia;)i+f^] are 
the smallest orthotopes aligned with the coordinate axes that contain AB and A~^B' , respectively. 
Since the above two sets of inequalities hold if and only if B and B intersect B' and B' respectively, 
the ^ direction is easily checked. In order to prove the <^= direction we make use of the fact that 
two disjoint convex polytopes can always be separated by the hyperplane supported by one d — 1 
dimensional face of one polytope. Thus, if AB and B' are disjoint, we can choose a face of B' or a 
face of AB. This respectively implies that BOB' = or S n S' = 0, and hence ends the proof. □ 



3.5 A priori estimates for a fully discrete LTP scheme 

In this section we establish a priori estimates for the fully discrete particle scheme consisting of 
an initialization step using the B-spline quasi- interpolation (3.3), and of a series of transport steps 
using the LTP transport operator ( 3.27 )-( 3.29), with no remappings. Thus, we consider 



and 



rn 
Jh 



n+1 



TJifJi for 



(3.40) 



with At = T /N . L°° convergence of this scheme will be established in Theorem 3.10 together with 
uniform bounds for the particle overlapping. A local (single-particle) error estimate will be shown 
as well, that will be of practical use in the adaptive multilevel LTP scheme described in Section |4j 
To express our estimates in terms of the smoothness of the velocity field u, we begin with local 
bounds for the characteristic flow. 

Lemma 3.6. Given a domain uj C M.'^, an integer m and two instants s,t G [0,t], we denote 

|u|m,(s.i.cu) sup |u(t', Olm.F^ 

where Fs,t' is the flow of u between s and t' , as above. Then we have 

\Fs,t - I\i,u, < ci,„exp(ci^„) and \Fs^t\2,u: < C2,« exp(ci,„)(l + Ci,„ exp(ci^„)) 
where c™^„ = 1* - s||u|m,(s,t,c^) form = 1,2. 

Proof. Rewriting (2.2) as dtFs,t{x) = u(t, Fs^t{x)), we obtain that for i,l — 1, . . . ,d, 

dtdi{Fs,t-iUx) =didt{F,M^) 

^ELidi'U,{t,F,,t{x))di{F,,t)i>{x) 

= T.t=idvu,{t,Fs,t{,x))di{Fs,t- I)i'{x) ^ diu,{t,F,,t{x)). 
In particular, using Fs^s = / we find 

d 



di[F,^t - lUx) = 



[Y^^i'^^^t' ^Fs,v^^))di{Fs^t' - I)v{x) + dMt\F,,t,{x)) 
i'=i 



At\ 



so that taking the supremum over x £ lo, summing over 1 = 1, 
I = 1, . . . , c? yields 



, d and taking the maximum over 



\Fs,t - /k^ < \u\us,t,u.) (I* - s| + - di') 

< 1^ - exp (|i - S\\u\x.(s,t.,u,)) 

where the second inequality follows from the Gronwall Lemma in integral form. This shows the 
inequality. In particular, using = 1 this shows that 



\Fs,tV,^ < 1 + |i - s\\u\xXs.t,L,) exp (|t - s||u|i,(s,t,t^))- 
Turning to the second derivatives, using again Fg^a = / we write 

d,MFsM^) = /.* [Eti-i5ii5,^ii»(i',^^M'(^))ai,(F,,,0;;(:i;)ai,(F,,,0;i(a^) 
+ Y.i'^idvu,{t!,F,^t,{x))diMF,^r)v{x)\ At' 



(3.41) 
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so that taking the supremum over x € uj, summing over li,l2 = 1, ■ ■ ■ ,d and taking the maximum 
over i = 1, . . . ,d now gives 



IF. 



t 



dt' 



< \u\2^is.t.u)iIjF,^t'\l^dt')exp{\t-s\ 

< C2,„(l + ci,„exp(ci^u))^exp(ci,„) 



u 



where we have used again the Gronwall Lemma, and (3.41 1. 
Corollary 3.7. For a given domain ui C M*^, we denote 

:= |ukcexp(At|w|i,5) and v^^ |u|2,5 exp(At|u|i,5)(l + Atz/i"_„), 
where ^ — (t" , t"^^ , lo) . Then we have 

In particular, denoting vf :— ^'"jjd yields 

ira^'ii < 1 + AiK\ ira^'b < Ati.^". 



□ 



(3.42) 



(3.43) 



Proof. The first estimate foUows from Fti^+i t{F^^{x)) = Ff^^tix) and Lemma 3.6 
foUows from the fact that F^" is a diffeomorphism. 



The second one 

□ 



Assumption 3.8. In view of Lemma 3.4 and Corollary 3.1 . one can expect that the error resulting 
from the linearization of the flow around x^'^^ behaves like 

\\J]:ix ~ xl+') ~ {{F:1)-\x) - 4!) II < Atihpfil + At< + h-'el,f'-'\,.^^ + 



so that it seems natural to ask that h 
such that 



-At 



At. In the sequel we will assume that h and At are 



h-'^e\t < aAt, 



(3.44) 



where a > is a fixed constant. Note that if an r-th order ODE solver is used to compute 
the numerical flow, the above condition reads At < Cah^ with a constant C depending on the 
smoothness of u, typically through |u|r_|_]^ (jn jn+i jfd-) . 

For the subsequent analysis it will be convenient to introduce the following measures of the 
velocity smoothness, 



K 



"1, 

n 

2,{h,k) 
n 

'3,{h,k) 
n 

'4,(h,k) 



I n 



hda, 
+ da, 



- 2d{l + Atn'l f^f^ i^-^y "'2,(/i,fc)' 
= 2d2(l + hAtnl^,^,^)\l + ^tnl^n,k)?('^-'^^liKk) 



(3.45) 



as well as space-global versions k"^ :— Vi +hda, . . . , and finally the time-global counterparts k^.^, 
z = 1, . . . , 4, obtained by replacing At with r, and using 



I'l = |m|i.j exp (Af|u|i^^) and |w|2,^ exp (At|u|i^^) (1 -I- Aii^i) 



(3.46) 



where ^ = (0, r, M''). Equipped with the local measures ( 3.45[ ) we derive the following estimates 
from Lemma 13.41 



Corollary 3.9. Provided h and At satisfy (3.44) and the additional mild condition 



the finite difference scheme (3.28) computes an invertible matrix JJ} satisfying 

\det{J^)\-' <{9^^)-' <l + hAtKl^^^,^, 



(3.47) 



(3.48) 
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therefore the operator is well defined by (3.27 1 -(3.29 1. Moreover, the linearized flow based on 
the approximated backward Jacobian JJ} satisfies the local estimate 



sup 



\\Jki^ - ^T^) - {ip:.r\^) - ^fc)iioo < h'^tciM 



(3.49) 



for all p>0, with Clk{p) = P<(/«,fc) + ^^'l(h.k) + 5/''<(Fji)-i(B„(."+\p/.))- 



Proof. According to (3.42), we have 



and 



(3.50) 



For /lAtKg < 1, we then infer from (3.30) that 

\d-l 



det( J,") >ei>l- hd{l + Atnl^,,^,^r-'Atnl^,^^,^ > 1 - |At<(,^,) > |Ai<(,_,) 



which shows the invertibihty of JJ} and also gives (3.481, hence the first claim. The local estimate 
(3.49) follows then from Lemma 3.4 and the definitions (3.45). □ 

We are now in position to state a priori estimates for the fully discrete LTP scheme. Again, 
we emphasize that although the following result is established for B-spline particles, the same 
arguments apply to more general approximation settings such as continuous finite element basis 
functions, see Remark 13.31 



Theorem 3.10. Let h > be lower than a fixed but arbitrary h, and At :— t/N be such 
that conditions (3.47) and (3.44) hold for some fixed a > Q. Then if the velocity field u is in 
L°°([0, r]; W'^'°°{W'-)), the numerical densities computed by the particle scheme (3.40) converge 
towards the solution f of equation (2.1), with 



IIA"-/(i")llL~ <Mct||/°||l^+ca|/°|i) for n = 0,...,7V, 



(3.51) 



with constants ct, ca independent of h and At, that are specified in the proof. The particles 
'^'hk~ ^h{D]!:{- — x^)) composing have uniformly bounded supports, 



^^PP{fh,k)^Bao{x'k,hph) with Ph = Cpexp{T{!^K3^h + Ki^h)), 



(3.52) 



in particular, ph — > Cpexp(Ti^i) as h ^ Q. Moreover, the overlapping constant of the particles 
Qh ■■= sup„<^^^gRd #({fc € Z"^ : ifiUx) ^ 0}) satisfies 



0/i < (2exp(Ti/i)(/9/i + hTa)Y {2cp)'^ exp{2dTi'i) as /i — ^ 0. 
Finally, the single-particle transport errors can be estimated by local indicators, 

{T]i-T-^WU < h^-''At\^\i\m\o.Cl,{pl+}) 
and where 



(3.53) 
(3.54) 



where CJ^ is as in Corollary 3.9 



,hAt + cpll (i?^)-i|U max { (1 + hAtKl^f^^,^)^! + At<^(^^,)), 1 + AK,e^ , } • (3.55) 
Proof. Following the main lines of the error analysis in Section |3.2[ for all n and k we let 

^n+l . /T^n rjin \,„n ,n-\-l rpn , „n 

■= - T^^Wh^k = fh,k ~ T^y,fh,k 
denote the particle transport error, and we define the corresponding overlapping constant as 

e„:= sup #{{keZ'' ■.el^,{x)^0}). (3.56) 



Thus, using the bound |w/c(/°)| < ^''llall^i satisfied by the weights in (3.3) we have 

\\{tj: - Tc" )/h lU- = II E ^^^fytkWL- < h^MU^wf^u^e, sup \K+}h-'. (3.57) 
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Next we claim that the local transport errors have their supports bounded with 



(3.58) 



where f. is defined as in (3.55). Using Corollary 3.9 this gives 

"h,/c = sup 

< \^h\i\m\oo 



sup 



< h^-'At\^U\\D-\\^Cl,{p-+}) 



where we have used \(ph\i 
estimate we observe that 



= h ^ '^\'p\i- This is our local estimate (3.54). To derive a global 
\{K-T:^)fH\\L^ < Ath(3{h,At)\\fU^ 



holds with 



sup 



{\m\ooCi,{pix')}. 



(3.59) 



It follows that the global error f^+^ - /(t"+i) = T^\/,7 - T^" /(i") is bounded by 



ii/r'-/(i"+^)ii 



<h{f3{h,At)\\fh^+CA\fU) 



(3.60) 



where we have used the fact that is an exact transport operator in the second inequality, and 
estimate (3.5) with g = 1 in the last one. In particular we can take 

CA = ^{mp + Cp)^{{2cp)'''\\ip\\L'^\\a\\ii + 1). 

It order to bound /3{h, At) independently of h and At we next observe that the deformation 
matrices involved in the scheme (3.27) verify 



with 



det(J™)rf(J™) ^ in the conservative case, or 
(J/T)^^ the non-conservative case. 



Using estimates ( |3.48[ ), (3.34) and (3.501 we then observe that 



|det(J," 



9^)-' < (1 + hAtK^^^,^,^)-^, WJnio. < t^T,k < (1 + ^t<,iH,k))- (3.61) 



It follows that 

||(i?fc)"'||oo < il + hAtK3,h)^il + AtniMT<cMr{'j^i3M + ni,h)) (3.62) 

holds for n < TV, in both the conservative and the non-conservative cases. To obtain a bound for 
|l(iI'j!)~^||oo we cannot repeat the above steps because our estimate (3.37) for ||(J™)~^||oo is too 
weak. However we can use the cofactor formula on directly. Indeed, by construction we have 
det(£'j?) = 1 in the conservative case, and in the non-conservative case we infer from (3.481 that 

|det(-DI?)| < n |det(jr)r' < (1 + hAtn^.j^r < expirhn^^h) 

m— 



and we observe that (3.62 1 becomes || (f^) ^||oo < exp{TKi^h)- Using the cofactor formula we then 
find in both cases that 



PI?IU < ddetiDmDkr'W'l^' < dexp(r(/i^3,h + Aii,/.(rf- 1))). 



(3.63) 
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This allows to bound the particles support. We indeed observe that 



El, supp«,) - + {D^r\B^{0,hcp)) C B^{xl,hcp\\{D]ir^^), 



(3.64) 



which shows the uniform estimate (3.52). From (3.64) we also derive that 
which allows to localize the particle transport error e^^^ with 



Using (3.44), (3.61) and ( 3.42^ this proves our claim (3.58). Turning next to the particle overlapping 
we consider an arbitrary x G M.'^, and for a given n < N we denote a;" = x and recursively set 
a;™ := (F™)-i(a;'"+i) for m = n - 1, . . . , 0. Using (|3l3| and ([Hi]) we see that 



\\x°-xl\\^<\iF^^)-M'x'-F^^{xl)\\^ 

< {l + Atu,){\\x^-xl\\^+elt) 

< . . . 

< {l + AtiyiY'iWx-x'^Wc^+nh^Ata). 

In particular, for k such that ^.(a:) ^ we have ||x — a;^'||oo < hph according to (3.52), hence 

fc||oo = - a;°||oo < exp(ri/i)(p/i + hra). 



\h-'x° 



This proves the estimate ( 3.53[ ) on the particles overlapping. As for the particle transport errors 
(3.56) we have 

< 29^, 



given that the particle overlapping is not increased by the exact transport operator. From (3.62) 
and vi < Ki^h we then observe that 

Pl,k < (^hAt + ph< aJiT + pj^ =: p-f^ 

holds for all At <t, n < r/At, fc e and h <h. In particular, the above discussion gives 

Pih, At) < Way |(^|i2(2exp(ri/i)(p;,))''dexp (t(/ik3_,- + K^-^^{d - l)))Cf,{pr,), 

where we have set Ch{p) = pi^i,h + i^3,h + \p^^2 - This latter bound provides us with a constant 
Ct for the global estimate (3.51 ), which completes the proof. □ 



4 Adaptive transport with multilevel particles 

To illustrate the flexibility of our approach we now describe one adaptive multilevel version of 
our LTP scheme, where the goal is to save computational effort without deteriorating the order 
of accuracy of the simulations. To implement that principle we use a hierarchy of particles with 
dyadic scales h = 2~^ , corresponding to integer levels j = jo, . . . , jmax- For notational simplicity we 
will label with j the objects that were previously labelled with their resolution /i, or that depended 
implicitly on it. For instance, particles with resolution h = 2"^ will now be denoted 

= )), 



instead of (^^ ^ as in equation (2.9 1. 

In Section |4.1| we first present one B-spline version of the common hierarchical approach for 
building adaptive approximations with prescribed error tolerance e, and in Section |4.2| we sug- 
gest a local filter that makes the resulting adaptive approximations positivity-preserving, without 



noticeable loss of accuracy. In Section 4.3 we then describe a dynamic strategy for refining the 
particles in the course of their transport, so that the associated transport error is on the order of 
the prescribed tolerance. 
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4.1 Adaptive approximations with multilevel B-splines 

To build adaptive B-spline approximations of a given function g, we reformulate the quasi-interpolation 
operator (|3.3| in a hierarchical setting. Specifically, writing 



A, : g ^ J2 "^jAaWlk with w.^kig) ■■= 2'^^ ^ a, 5(2"^ (fc + 0), (4.1) 

we replace the approximation g « Aj^_^^g by a recursive process where we start from gj^ := Aj^g 
and define gj by correcting for j = jo + 1, • ■ • , jmax, i.e., 

In such a framework, adaptivity is classically obtained by restricting the corrections in the regions 
where the residual is larger than the prescribed tolerance. This idea can be realized by discarding 
small particles: indeed, using 

Ay.g^ ^ w,A9)^lk with Q{e):=2-''^e{{2cpfML^r' 

'=:|"'3,fc(s)l>Cj(e) 

gives an acceptable error, as 

\\{A^ - A,)g\\L- = II E "^^^9)^14^^ < 2*(2c,)''||^|U^C,(e) - e. 

k:\w,_^{g)\<C,{e) 

Moreover, the adaptive correction process only needs to be applied in the regions where the residual 
is significant. Specifically, let 

n+ {x : |(.g - gj-i)ix)\ > c+e} and K+ {k : + mp[-l, 1]'') n n+ ^ 0} 

where c'^ := ((2cp)''||(y9||L°° ||a||fi)^"'". It is then easily seen that 

- 9]-i)\ < 0(e) for k ^ K+ . 

In practice the set may not be readily computable, but it is possible to find some reasonable 
approximations K* that contains it. The resulting scheme for adaptive quasi-interpolation with 
multilevel B-splines reads then as follows. 

Algorithm 4.1 (Adaptive B-spline quasi-interpolation, A^ : g i-^ .9imax)- 

1. Initialize the algorithm with gj^-i :— and pre- activate every coarse particle (i.e., make it a 
can 



didate for a possible selection) by letting K*^^ := Z . 
2. For j = jo, ■ • • , jmax, do: 



(a) for k £ K* , compute w* ^. = Wj,k{g — gj-i) as in (4.1|, then set 

9j — 9j-i + E ^Ik^lk where := {k e K* : > Cj(e)}; 

(b) ^f i < jmax, pre- activate particles at the finer level for possible selection by determining 
Kj^i (see Algorithm below); 

(c) if the option is set, apply the positive correction filter (see Algorithm 4-4 below). 

To determine the set Kj+i of candidate particles, a naive approach would evaluate local values 
of the residual g — gj. However this is somehow redundant with the computation of the weights 
f^j+i k '^^^t level, therefore it should be avoided. Instead, we observe that Kj+i can be built 

as follows. 

Algorithm 4.2 (Pre-activation of particles in Kjj^i, for j = jo, ■ • • , jmax — !)• 
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1. Initialize the set with K^j^i = 0- 

2. For k e K* , do: 

(a) if k Kj has been selected, add {2k + I : \\l\\oo ^ '^ip + 2cp} to K*^^ 

(b) otherwise, evaluate the residual rj^k '■= \\9 ^ 9j-i\\L°°(B^(2-ik,2-i-'^)) '^'^'^ <^dd {2fc + / 
II^IU < rnp} to K*^^ if rjM > c+£. 



Proposition 4.3. The sets built in Algorithm \4-2\ do satisfy 

Kj + l C — JOi • ■ • I Jmax — 1- 

Proof. Take k G K^^-^. By definition there is an x sucli tliat \{g—gj){x)\ > c+e and ||2-'+^a; — /c| 



TOp. Then if k has not been added to K*_^_i on step 2a it is easily seen that \\2^x — fc'||oo > Cp holds 
for all k' G Kj. In particular we have \{g — gj-i){x)\ — \{g — 9j){x)\ > c^e, hence x G fi^. Take 
now k" such that \\k" — 2^a;||oo < | (this is possible up to choosing another x close to the first one, 
as g — 9j-i is continuous). This gives rj^k" > c^s and k" G Kj' , by definition of the latter. Now, 
it also gives \\2k" — fc||oo < 2||fc" — 2-'x||oo + ||2^+-'^a; — fc||oo < riip + 1. Since we can assume as an 
induction hypothesis that K^ C K* (it is valid for K*^^ = Z''), k must have been added to Kj_^_^ 
on step[2b| □ 



Finally we note that the main loop in Algorithm 4.2 can be merged with step 2a from Algo- 
rithm 4.1 In particular, the local residuals rj^k can be evaluated while computing the weights 



w 



4.2 Positivity preserving hierarchical approximations 

If the target function g is nonnegative, one may want its particle approximation to be nonnegative 
as well. However, it is easily seen that the above scheme is likely to yield some negative weights. A 



first reason for this is that the single-level scheme (4.1 ) is not positive for p > 1. A second reason 
is due to the hierarchical framework: wherever gj-i is below the target, the residual takes negative 
values which are likely to be approximated by negative level-j particles. 

The latter issue can be addressed by locally decreasing gj_i in the neighborhood of a negative 
weight Wj^k, so as to increase the residual and hopefully correct the negativity of Wj^k- Since 
B-spline particles satisfy a refinement equation 

'Pj^k' = '^lf3' + l-2k'+l, (4.2) 

l|;|| = <Cp 

it is possible to do so by refining coarse particles which contributions to ipj^k will increase its weight. 
An attractive feature of such corrections is that they do not require new evaluations of the target 
function, or iterative approximations of updated residuals. Let us specify one such algorithm. 
Because it may be necessary to refine particles over several levels, it is useful to derive multilevel 



scaling relations from (4.2 1, 

Vj',k' = '^i^^y^3'+s,2'^k'+i (4-3) 

ll'lloc<sa 

with multilevel refinement coefficients satisfying 



a 



and initialized with ad) = a, s, — Cp. Indeed, we can observe that if all the particles at the level 
j — S, S > 1, were to be refined up to level j, the resulting increase to Wj^k would be 



fc'ez 



19 



Therefore, letting Sj'{j, k) := Wj^k + k) we define the correction level associated to 

some given negative weight Wj^k as 

j_ (j, k):=l < ^ ■ ^^-^^ - (4.4) 

' 1 jo otherwise. 

It corresponds to the coarser level from which neighboring particles should be refined in order to 
correct the weight Wj^k- Obviously, it is not necessary to refine every such coarse particle, since 
only those satisfying ||fc — 2^k'\\ao < ss do contribute to (j, k). Moreover, we do not need to fully 
refine them up to level j. Instead, the same increase to Wj^k can be obtained by refining them 
one level at a time, and only considering those contributing to (j, k) in the process, namely those 
{j',k') with \\2^^^ k' — k\\oo < Sj-^j'. Thus we see that these corrections are indeed local. 

We also observe that for p > 1, it may happen that some weights remain negative despite the 
corrections. In this case we suggest to simply discard those weights after correcting them. Clearly 
this may deteriorate the approximation accuracy, but in a hierarchical framework one can hope 
that finer layers of details will essentially correct the resulting errors. In practice indeed we have 
observed such a behavior, see Section[5.2[ Let us now summarize the above correction filter, to be 



applied to the weights iWj'.fc, j' = jo, ■ ■ ■ ,j, computed on step 2a of Algorithm 4.1 Note that here 
we compute the correction levels before refining the contributing particles so as not to break the 
possible symmetries in the particle grid. 

Algorithm 4.4 (Positive correction filter). Let KJ~'^ := {k e Kj : Wj^k < 0} denote the indices of 
the active particles to correct at level j . Then, 

i- if i > jo, 

(a) set jcoi{j, k) < j as in (4.4), for all k G KJ''^ ; 



(b) then, refine the contributing particles. Namely, for k G K'^^ and j' — jcor(j, k), . . . ,j — \, 
refine every particle (j',k') such that \\2^^^ k' — fc||oo < ^j-j' setting 

Wj'+i,2k'+i ■■= w-j'+i,2k'+i +a-iWj>^k' for all \\l\\oo < si, and wj'^k' ■= ; 

2. for all k G K^'^ , set wj^k '■— max{0, Wj^fc}. 

Remark 4.5. We observe that there is no need to pre-activate new finer particles after running the 



above correction algorithm, i.e., the set K*^-^ built in Algorithm 4.2 still verifies Proposition 



4.3 

Indeed, the refinements on step lb do not change the approximation gj, nor the residual, i'he 
only change to gj occurs when discarding a negative weight Wj^k- There, a significant error may 
be introduced but since k is in Kj, a finer patch has already been pre-activated on step 2a of 



Algorithm 4.2 



4.3 Dynamic particle refinement without remapping 

In order to achieve some prescribed accuracy when transporting the multilevel particles 



j=jo fceZ'' 

it is in general necessary to refine some of them over time. Indeed, the errors induced by the discrete 
transport operator depend on the local smoothness of the flow, and on the level of the particles. 
However, at the initialization step their resolution is essentially driven by the local smoothness of 
Rather than follow a conservative approach and automatically refine patches of particles to 
resolve the emerging features that may appear in the solutions, we propose to use the local error 



estimate (3.541 to determine which particles are admissible for discrete transport, and which need 
to be refined. Given that both the number of levels and the number of overlapping particles per 
level are uniformly bounded, we have \\{T" ~ ^ox)/jL^^ IU~ ^ ^^Pj,k\\''^j',k^]t^\\L°° hence a particle 
w'^kV^j k ^^i*^ admissible for transport if it satisfies 

\wlk\v7,k<CTe, (4.5) 
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where ry"^, is a computable estimate for the transport error, see (4.7) below, and Ct is an ad-hoc 
constant. 

We next observe that, although it is possible from (4.2) to represent exactly a deformed particle 
by a few finer ones, using 



ii;|ioo<cp 



)) 



where 



2 ^ ^{D'^k) refining inadmissible particles in such a way can result in a 



significant increase of their total number. Assume indeed that every level-j particle needs to be 
refined in /". Since there is no reason why we should have D^^ — fc' ^'^^ fc / ~ fc' v f*-"^ ^' 
and I' such that 2k + I — 2k' + I' , replacing every ip^j^. by its exact finer representation is likely to 
add ~ 2'*' (2cp + l)'* new level-(j + 1) particles, which is much more than the particles involved by 
a uniform discretization at level j + 1, and similar phenomena may happen on further time steps 
as well. 

For that reason we propose to refine inadmissible particles v?"^, in a "retroactive" fashion. 
Namely, we represent them by the level- (j -1-1) particles that would have been present in /" if 
the original (structured) (p'j had been refined from the start, that is, at the last remapping step. 
Note that some of those fine particles may already be present in /" , and in such a case it suffices 
to update their weight. In particular, this approach allows to refer to any particle by its multilevel 
space-time index without ambiguity: for all j, k,n we indeed have = ¥'j(-D"fc(- — a;"fc)) with 
particle center and deformation matrix given by 



n-l 

. . . ^«o)(2-.fc) and Dl, = Dj^^'j;^' ^ (4.6) 



where no is the last remapping step and J™^ is defined by ( 3.28 )-( 3.29). 

It remains to give one error indicator Ty"^. for the local admissibility condition (4.5). In the 
numerical tests shown in Section [5] we have used local estimates 



d 

max V|A>,(r)(x^,)|, \u\l^„ 
1=1 



max y 

l<i<d 

/l,i2 = l 



Af^,^u,(t")(x;^,)l, 



based on the finite differences 

A/w(t")(x) = 2\u(e,x^ 2-^-^ei) - u{t",x- 2-^-^ei)) and Af^ = 
We mimic the local smoothness measures from Lemma [3^ and Corollary |3.7| with 



Af A' 



i'um ■■= l"l2..^,,exp(At|u|t^,„J(l + Atz>J^(^.,)) 



and define coefficients /t"j^^ i 



1, . . . , 4, consistent with (3.45). Following (3.54| we finally take 

:= 2^(''-i)A%|ip;^,||ooq:fe(p;^r) 

1 ^2,- 



(4.7) 



where we have set C^^^ip) — P^lij.k) + '^lij.k) + iP Kum) ^^'^ 



Note that neglecting the small terms in (4.7) gives 

r,l, « V^^~'^At\cp\,\\Dl,\\^ (ip2 + + 2<fa{dp + 1)' 



where p = Cp|| (1?!^^^)- 



/5"J and a is the constant from assumption (3.44). We summarize 



our adaptive transport algorithm for multilevel particles as follows. 
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Algorithm 4.6 {T^ : ^ f^+l). For j = jo, ■ ■ 
do: 



and then for k G Z'^ such that ^. 7^ 0, 



*/ j = imax or if the admissibility condition (4.51 is met, then the particle is transported by 
setting wj+^ := wl^, x'Jj:^ := F"(a;J J and DJ+^ := D^i^J^,^ according to ^^-^^ ; 



2. otherwise, it is dynamically refined by adding (Tiw^j ^. to w^j^-^ 2k+i f'^^ every I such that \\l 



< 



Cp (which may involve the activation of new particles according to (4.6)^ and by finally setting 



Remark 4.7. One disadvantage of tlie above strategy is the need to apply the numerical flows from 
past time steps, namely F"°, . . . , F"~^, when refining a particle at time step n. However it is 
useful from a theoretical point of view. Practical procedures involving local approximations of the 
time-integrated flow are currently being tested and will be presented in a forthcoming article. 

5 Numerical experiments 

In this section we compare the numerical performances of the different particle methods, using 
Leveque's swirling velocity field [21j in two dimensions. 



UQ{t,x) = — sin^(7ra;o) sin(27ra;i)(7(i), ui{t,x) — sin^(7ra;i) sin(27rxo)g(i) 



(5.1) 



where g{t) — cos(7ri/T), t e [0, T]. The corresponding flow is easily pictured from the space pattern 
shown in Figure [Tj and from the symmetry of g with respect to T/2 which reverts the solutions to 
their initial state at t = T. We will take T = 2.5 which corresponds to a moderate stretching at 
t = T/2 (see Figures [I] and |6] below) and a fixed time step At = T/100. 

Although the time symmetry may simplify the measure of the numerical errors, we note that 
by considering only the final accuracy one misses the intermediate errors resulting from the inac- 
curate transport of the particle shapes, which yields biased estimates for particle methods with no 
remappings. Therefore, we shall only use the time symmetry in Section [5.2| when assessing the per- 
formances of the uniform and adaptive versions of the linearly transformed particle (LTP) scheme, 
with at least one remapping. In Section 5.1 where the traditional smoothed particle method (TSP) 
is compared to the forward semi-lagrangian (FSL) and the LTP schemes, we will only solve the 
equation on the half interval [0,T/2], and hence measure the numerical errors at the instant of 
maximal stretching, using a reference numerical solution obtained with significantly higher space 
and time resolution. 
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Figure 1: Leveque's swirling velocity field (5.1 1 shown in the computational domain [0, 1]^ at time 
t = 0. 
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5.1 Numerical comparison of the TSP, FSL and LTP schemes 

To compare the LTP scheme with the classical TSP and FSL methods reviewed in Section [2] we 
consider the two following cases. 

• A smooth hump centered on c = (0.5, 0.75) given by 

/"(x) = i(l+crf(a + &||a;-c||2)) with a = 3.43, 6 = 2L43. (5.2) 

Here ||-||2 is the Euclidean norm and erf(s) = ^ e^^ dy is the standard "error function" 
that smoothly spans [—1, 1], see Figure [2j 

• A discontinuous Zalesak's slotted disk centered on c' = (0.5,0.7), 

f(x) = iJ(0.15 - llx - c'||2)(l - (1 - H{\xo - 0.5| - 0.02))(1 - H{xi - 0.8))) (5.3) 
where H{s) ~ X{s>o} denotes the Heaviside step function, see Figure [i] 

In Figures [3] and [5] the relative errors are plotted versus the average number of active particles 
a.t t — T/2, in order to avoid a biased measure for the TSP method as explained above. For 
the smooth hump case we measure the errors in and for the discontinuous Zalesak disk we 
use the norm. Here all three methods have been run with B-spline particles of degree p = 1 
or 3, initialized (and remapped, in the FSL and LTP cases) with the quasi-interpolation scheme 



(3.3 1 of corresponding order. For the TSP method we show the effect of varying the overlapping 
exponent q such that e = /i*, see Section [2T] and for the FSL and LTP schemes we plot the runs 
corresponding to decreasing remapping frequencies, i.e., increasing values of the number TV,- of time 
steps between two remappings. Observe that the higher q or N^-, the cheaper the simulations. In 
every set of curves we have used a thick line to emphasize the cheapest run: for the TSP scheme 
that corresponds to the case q = 1 (particles radii are proportional to the initial meshsize) and for 
the FSL and LTP schemes that corresponds to the case with no remappings. Note that the FSL 
and TSP method coincide with such parameters. 

From these plots we make the following observations. First, we see that a necessary condition 
for the TSP method to converge is indeed that particles present an extended overlapping. That 
is, the ratio e/h must go to +cxd as h goes to 0. Moreover, the convergence is relatively slow 
and is not much improved with third-order B-splines, although they help for moderate particle 
overlapping {q 0.8). This is not surprising since the theoretical convergence analysis requires 
vanishing moments for high-order accuracy, which the B-splines do not have. 

Second, our numerical tests confirm the announced practical condition for the FSL scheme to 
converge: remapping frequency must be high. On the plus side, we note that on the smooth test 



case (5.2), the FSL scheme exhibits significantly higher convergence rates than the TSP method 
when particles are used. However, to our view the method is severely hampered by the following 
dilemma: for values of A^r close to 1, problems with sharp edges such as the discontinuous Zalesak 
disk give rise to a strong numerical diffusion that can be seen (especially with B3 particles) from the 
increase of the number of active particles and from the deterioration of the accuracy ; decreasing the 
remapping frequencies helps reducing that effect, but doing so always leads to a loss of convergence 
(the lower the frequency, the sooner the stagnation). 

Finally, we see that the LTP scheme always converges as expected from our analysis, including 
in the cheap, non-remapped runs, although remappings do improve the accuracy in the smooth 
test case (Figure |3]). Moreover, the convergence rates are significantly higher than with the TSP 
method, and in the smooth case the benefit of using ^3 splines is obvious, just as in the FSL 
method. In the discontinuous case (Figure [5| we observe that lowering the remapping frequency 
strongly reduces the numerical diffusion (especially with particles), and in fact also improves 
the numerical accuracy. The striking result is that the loss of convergence observed in every FSL 
run using a low remapping frequency is completely suppressed in the LTP runs. 

5.2 Numerical study of an adaptive LTP scheme 

To assess the ability of our adaptive particle strategy to improve the computational efficiency, we 
compare its convergence curves to those of the uniform LTP scheme, using the two following cases. 
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Figure 2: Profile of the exact solution at t = or T (left) and r/2 (right) for the smooth hump 
(pi. 
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Figure 3: Convergence curves (relative errors at i = T/2 vs. average number of active particles) 
for the smooth hump test case shown in Figure [2j using the TSP, FSL and LTP schemes with B- 
spline particles of degree 1 and 3. Triangles represent slopes of 0.5 and 1. 
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Figure 4: Profile of tlic exact solution at i = or T (left) and r/2 (right) for the Zalesak's slotted 
disk (lO). 
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Figure 5: Convergence curves (relative errors at t — T/2 vs. average number of active particles) 
for the Zalesak's test case shown in Figure [4j using the TSP, FSL and LTP schemes with B-spline 
particles of degree 1 and 3. Triangles represent slopes of 0.5 and 1. 
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• The smooth hump (5.2) shown in Figure [2] 



A smoothed version of the discontinuous Zalesak's slotted disk, designed so as to get solutions 
of highly non uniform smoothness that still can be approximated in the uniform norm. 
Indeed our adaptive refinement strategies are aimed at balancing the local errors in this 
norm and they have no reason to exhibit better convergence properties when the error is 
measured using another norm. The initial data /° is then obtained by substituting the 



discontinuous Heaviside step function H(s) — X{s>o} (5.31 with a smooth approximation 
H^{s) — ^{1 + erf(s/e)), where the transition zone has approximate diameter of 4e. In the 
numerical tests we take e = 0.01, see Figure [6j 

In Figure [7] we plot for both cases the L°° convergence curves obtained with the uniform LTP 
scheme using various values of TV^ as above, and those obtained with our adaptive particle method 
using the optimal value TV]. = 25 as found from the uniform runs, and levels j — jo, - ■ ■ , Jmax with 
jo — 5 and jmax = 9. Since we now consider particle methods with at least one remapping step, 
numerical accuracy is measured on the final time t = T where the exact solutions revert to their 
initial state. The prescribed tolerance e varies in order to get different points in the "adaptive" 



convergence curves and we have fixed the constant Ct in Condition (4.5) to 2, to optimize the 
numerical performances. However we note that a proper study of how the parameters Ct and N,. 
should be set in the general case remains to be done. 
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Figure 6: Profile of the exact solution at t 
Zalesak's slotted disk. 
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Figure 7: L°° convergence curves for the smooth hump (left) and the smoothed Zalesak's slotted 
disk (right) in Leveque's swirling flow. The relative L°° errors are plotted vs. the average number 
of particles. Here the adapt curves correspond to adaptive runs and are obtained by varying the 
tolerance e. The other curves are obtained by varying the uniform resolution h, as in Sect ion |5.1[ 
In the right panel the runs labelled adapt + use the positive correction filter (Algorithm 4.4) in 
each adaptive remapping. Triangles represent slopes of 0.5 and 1. 
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particle centers at t=T/2, before remapping 



i. 
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particle levels at t=T/2, before remapping 



particle centers att=T/2, after remapping 
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particle levels att=T/2, after remapping 



Figure 8: Particle centers (left) and level maps (right) a,t t = T/2 for the smooth hump test case 
shown in Figure [2j before (top) and after (bottom) being remapped on a hierarchy of cartesian 
grids with resolution levels j = 5, . . . , 9 at i = T/2. Here the prescribed tolerance e is set so as 
to get an average number of particles corresponding to a uniform run of level j — 7, see Figure [9] 
below. 
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Figure 9: Final error distributions obtained for the smooth hump test case with a uniform run of 
level j = 7 using about 10164 particles in average (left), and the adaptive run shown on Figure [s] 
using about 10498 particles in average (right). 
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Figure 10: Particle centers (left) and level maps (right) for the smoothed Zalesak's slotted disk 
shown in Figure [6j before (top) and after (bottom) being remapped on a hierarchy of cartesian 
grids with resolution levels j = 5, . . . ,9 at t — T/2. Again the prescribed tolerance e is set so as 
to get an average number of particles corresponding to a uniform run of level j = 7, see Figure |11| 
below. 



error at t=Twith particles of level j=8 



error att=Twith adaptive multilevel particles 



I 



Figure 11: Final error distributions obtained for the smooth Zalesak's slotted disk with a uniform 
run of level j = 7 using about 11864 particles in average (left), and the adaptive run shown on 
Figure 10 using about 11964 particles in average (right). 
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From Figure [7] we see that using adaptive particles with dynamic refinements yields a clear 
gain in terms of active particles number. Moreover, in the Zalesak case where the numerical 
solutions are likely to strongly oscillate in the vicinity of the sharp edges, we also plot a series of 
runs obtained with our multilevel correction Algorithm 4.4 that enforces positive weights, hence 
positive particles. The resulting observation is that the numerical efficiency does not seem to suffer 
much from this local (and somewhat crude) positivity-preserving filter. 

To illustrate the behavior of the adaptive particle method, we also show on Figures [8| and [l0| the 
distribution of active particle centers and associated level maps x i->- max {j : J2k<£Z''\''^^ k'P^ k\ 7^ O} 
at t = T/2 (before and after remapping the particles). 

Finally, one relevant feature of our adaptive particle scheme is shown on Figures H] and [11 



for the smooth hump and the smoothed Zalesak disk, respectively. There we have plotted the 
final error distributions obtained with uniform and adaptive runs using similar numbers of active 
particles (the adaptive runs correspond to those shown in Figuresjsjand 10 1. In both cases we verify 
that resorting to adaptive particles does not only reduce the maximum error, but also balances the 
error distribution over the computational domain, as it should be. 



6 Conclusion and further work 

We have introduced a formal class of particle methods for transport problems with polynomial 
deformations of arbitrary degree r, and established their L°° convergence with order r, based 
on the smoothness of the underlying velocity field. In the first order case we have presented a 
fully discrete scheme for the resulting LTP scheme, along with local and global error estimates. 
Our preliminary numerical results demonstrate the improved convergence rate of the LTP scheme 
compared to the traditional "smoothed" particle method (TSP) with extended overlapping, and 
the need of lower remapping frequencies compared to the FSL scheme, leading to lower numerical 
diffusion and computational cost. On a practical level, implementing the FSL scheme is made 
simple by the fact that it only involves pointwise evaluations of the forward flow, a building block 
in any particle scheme. We also note that pointwise evaluations of the deformed particle collections 



can be made computationally efficient by following the procedure described in Section 3.4 

In a near future we shall present extensions of the LTP scheme that deal with non-linear 
transport problems arising in, e.g., plasma physics [B], and we will further investigate the properties 
of our adaptive multilevel scheme. 
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